4. Wave Equation

We consider the wave equation as typical example of second order hyperbolic equation. We have a second order derivative in time, and second order derivatives in space, but of different signs. This completely changes the behaviour of the solution.

\[ \frac{\partial^2 u}{\partial t^2}(x,t) - \Delta u(x,t) = f(x,t) \qquad \forall \, x \in \Omega, \forall t \in (0,T) \]

Since it is second order in time, we need two initial conditions, for the function itself and for the time derivative:

\[\begin{split} \begin{array}{rcl} u(x,0) & = & u_0(x) \\ \frac{\partial u}{\partial t} (x,0) & = & v_0(x) \end{array} \end{split}\]

We use a variational formulation in space and obtain the second order ordinary differential equation (ODE) in time:

\[ M \ddot u(t) + A u(t) = f(t) \qquad \forall \, t \in (0,T) \]

with initial conditions

\[ u(0) = u_0 \qquad \text{and} \qquad \dot u(0) = v_0 \]

4.1. Newmark time-stepping method

We reduce the second order ODE to a first order system. For this, we introduce the velocity function \(v = \dot u\):

\[\begin{split} \begin{array}{rcl} \dot u & = & v \\ M \dot v & = & f - A u \end{array} \end{split}\]

Now we apply the trapezoidal rule for time-discretization:

\[\begin{split} \begin{array}{rcl} \frac{u^j - u^{j-1}}{\tau} & = & \frac{v^j + v^{j-1}}{2} \\ M \frac{v^j - v^{j-1}}{\tau} & = & \frac{f^j+f^{j-1} - A (u^j+u^{j-1})}{2} \end{array} \end{split}\]

Expressing \(u^j\) from the first equation, and inserting it into the second we obtain

\[ M \frac{v^j - v^{j-1}}{\tau} = \frac{f^j+f^{j-1} - A (u^{j-1}+\tfrac{\tau}{2}(v^j+v^{j-1})+u^{j-1})}{2}, \]

and shifting the unknowns to the left we obtain a linear system for the update of velocity:

\[ (M + \tfrac{\tau^2}{4} A) (v^{j}-v^{j-1}) = \tau \left( \tfrac{1}{2} (f^j + f^{j-1}) - A u^{j-1} - \tfrac{\tau}{2} A v^{j-1} \right) \]

Although by introducing the system of ODEs, the size of the linear system to solve did not increase.

from ngsolve import *
from ngsolve.webgui import Draw
from time import sleep

from netgen.occ import *
rect = MoveTo(-1,-1).Rectangle (2,2).Face()
rect.edges.Min(X).name="left"
circ = Circle ( (0.7, 0.0), 0.1).Face()
shape = rect-circ
mesh = Mesh (OCCGeometry( shape, dim=2).GenerateMesh(maxh=0.03)).Curve(3)
Draw (shape)
Draw (mesh);

tau = 0.002
tend = 4
u0 = exp(-400*( (x)**2 + (y)**2))
v0 = 0

fes = H1(mesh, order=3)
u,v = fes.TnT()
mform = u*v*dx
aform = grad(u)*grad(v)*dx

m = BilinearForm(mform).Assemble()
a = BilinearForm(aform).Assemble()
mstar = BilinearForm(mform+tau**2/4*aform).Assemble()
mstarinv = mstar.mat.Inverse(inverse="sparsecholesky")
f = LinearForm(fes).Assemble()

gfu = GridFunction(fes)
gfv = GridFunction(fes)
gfu.Set(u0)
gfv.Set(u0)

scene = Draw (gfu, deformation=True, euler_angles=(-60,0,-10))
sleep (5)
with TaskManager():
    for j in range(int(tend/tau)):
        gfu.vec.data += tau/2 * gfv.vec
        gfv.vec.data += tau * mstarinv * (f.vec - a.mat * gfu.vec)
        gfu.vec.data += tau/2 * gfv.vec

        if j%10 == 0:
            scene.Redraw()
scene.Redraw()